﻿using System;
using System.Numerics;

namespace HydrogenAtom4 {
    class HydrogenAtom {


        private static double LaguerrePolynomial(int n, int k, double x) {
            if (k < 0) {
                throw new ArgumentOutOfRangeException("k");
            }
            if (n < 0) {
                throw new ArgumentOutOfRangeException("n");
            }

            if (n == 0) {
                return 1.0;
            }

            double y = 0;

            for (int i=0; i <= n; ++i) {
                double t = Util.Factorial(n) / Util.Factorial(i);
                t *= Util.Factorial(k + n) / Util.Factorial(n - i) / Util.Factorial(k + i);
                t *= Math.Pow(-x, i);
                y += t;
            }
            return y / Util.Factorial(n);
        }

        private static RealPolynomialOneHalfCoeffPrecision BuildLegendrePolynomial(int n) {
            if (n < 0) {
                throw new ArgumentOutOfRangeException("n");
            }
            if (n == 0) {
                // 1
                return new RealPolynomialOneHalfCoeffPrecision(new double[] { 1.0 });
            }
            if (n == 1) {
                // x
                return new RealPolynomialOneHalfCoeffPrecision(new double[] { 0.0, 0.0, 1.0 });
            }

            return BuildLegendrePolynomial(n - 1).MulX().Mul((2.0 * n - 1.0) / n)
                .Sub(BuildLegendrePolynomial(n - 2).Mul((n - 1.0) / n));
        }

        class LegendrePolynomial {
            private readonly RealPolynomialOneHalfCoeffPrecision mCoeffs;
            private readonly int mL;
            private readonly int mM;
            public int L { get { return mL; } }
            public int M { get { return mM; } }
            public LegendrePolynomial(int l, int m, RealPolynomialOneHalfCoeffPrecision coeffs) {
                mL = l;
                mM = m;
                mCoeffs = coeffs;
            }

            public Complex ValueOf(double x) {
                return mCoeffs.ValueOf(x);
            }
        }

        LegendrePolynomial mLegendrePolynomialDifferential;

        public int L { get { return mLegendrePolynomialDifferential.L; } }
        public int M { get { return mLegendrePolynomialDifferential.M; } }

        private void BuildLegendrePolynomialDifferential(int l, int m) {
            if (Math.Abs(m) > l) {
                throw new ArgumentOutOfRangeException("l");
            }

            mLegendrePolynomialDifferential
                = new LegendrePolynomial(l, m,
                    BuildLegendrePolynomial(l).NthDifferential(Math.Abs(m)));
        }

        private Complex AssociatedLegendrePolynomial(double x) {
            int l = mLegendrePolynomialDifferential.L;
            int absM = Math.Abs(mLegendrePolynomialDifferential.M);
            Complex v = 1.0;
            if (absM != 0) {
                v = Math.Pow(-1.0, absM) * Util.PowComplex(1.0 - x * x, absM, 2);
            }
            v *= mLegendrePolynomialDifferential.ValueOf(x);
            return v;
        }

        private Complex SphericalHarmonics(double θ, double φ) {
            int l = mLegendrePolynomialDifferential.L;
            int m = mLegendrePolynomialDifferential.M;
            Complex v = Util.PowComplex(-1.0, m+Math.Abs(m), 2) *
                Math.Sqrt((2.0 * l + 1.0) * Util.Factorial(l - Math.Abs(m))
                          / (4.0 * Math.PI) / Util.Factorial(l + Math.Abs(m)))
                * AssociatedLegendrePolynomial(Math.Cos(θ));
            var eimφ = new Complex(Math.Cos(m * φ), Math.Sin(m * φ));
            return v * eimφ;
        }

        public HydrogenAtom(int l, int m) {
            BuildLegendrePolynomialDifferential(l, m);
        }

        public Complex WaveFunction(double a0, int n, double r, double θ, double φ) {
            if (n < 1) {
                throw new ArgumentOutOfRangeException("n");
            }

            int l = mLegendrePolynomialDifferential.L;
            int m = mLegendrePolynomialDifferential.M;
            double ρ = 2.0 * r / n / a0;

            double v = Math.Sqrt(Math.Pow(2.0 / n / a0, 3)
                * Util.Factorial(n - l - 1) / 2.0 / n / Util.Factorial(n + l))
                * Math.Pow(Math.E, -ρ * 0.5) * Math.Pow(ρ, l);
            return v * LaguerrePolynomial(n - l - 1, 2 * l + 1, ρ)
                * SphericalHarmonics(θ, φ);
        }


    }
}
